A first pass at a definition and implementation of a “nearly neutral” (maybe - this is the term we started with but it may not actually be applicable) model of the evolution of a parasite community within a single host population.

Start with a host population of size \(N\) and subpopulations \(\{n_i\}\) of different pathogen strains (the length of the vector will be dynamic …). Each strain has an associated \(R_{0,i}\) (the vector could be defined in sorted order - I don’t know if this would make either the mathematical description or the computational organization more convenient). For practical purposes, we will use the per-contact probability \(\beta\) instead, where \(R_0 = \beta N/\gamma\) and \(\gamma\) is the recovery rate (i.e., probably per time step), set equal for all parasites for the time being (although later we can imagine allowing this to mutate as well). Consider a discrete-time SIS model with mutation in the transmission probability.

  1. The first update step is infection: uninfected individuals (\(S=N-\sum n_i\)) remain uninfected with probability \[ \prod_i (1-\beta_i)^{n_i} \] (i.e., a Reed-Frost update). Those individuals that do become infected are subdivided with a per-strain probability of approximately \[ p_i = \frac{n_i R_{0,i}}{\sum_j n_j R_{0,j}} \qquad , \] i.e. a weighted fraction of prevalence. (Side question #1: what is the proper probability of being infected first by an individual of type \(i\) in this case? …) If we use the approximation, we could make it stochastic and discrete by drawing a multinomial probability … or, even more approximately, draw Poisson deviates with mean (\(S/N n_i R_{0,i}\)) for each infection type …
  2. In the second update step, all previously infected individuals recover with probability \(\gamma\) (might need to think carefully about the consequences of this particular update rule [infect first, then recover] to see if there are adjustments that need to be made to the relationship between \(\beta\), \(\gamma\), and \(R_0\) ….).
  3. Consider mutation of new infecteds. (Steps 2 and 3 can occur in arbitrary order.) Mutation occurs with per capita probability \(\mu\) (i.e., per new infection). If mutation occurs we pick a new \(R_0\) (equivalently, a new \(\beta\)) with one of the following rules:

\(\epsilon\) is drawn from a Normal distribution with mean \(a\) (\(<0\)) and standard deviation \(b\). (We may have to be careful about scaling: if we have discrete time scales then \(\beta\) may have to move/be defined on the log-hazard or logit scale instead.) 4. Initial conditions: somewhat arbitrary, but possibly a single strain with \(R_{0,1}=2\) and \(n_1=0.5\) (which is the equilibrium condition for a single strain).

Side question #2: what is the effective population size for this type of infection system? Logically, it seems that \(N_e\) should be 1 when \(I=1\) (because only one individual is contributing to the population) as well as when \(I=N-1\) (because there is only one susceptible left), and presumably maximum (\(=N/2???\)) when \(I=N/2\)

Preliminary results

Notes

These are the distributions of \(\beta\) and \(\log_{10}(R_0)\) (the latter certainly needs to be taken with a grain of salt …)

Latin hypercube runs

Run simulations across a broad range of parameters:

var min max
gamma 0.100 0.500
mu 0.001 0.100
mut_mean -4.000 -0.500
mut_sd 0.100 3.162
N 3.000 1000.000

For each run, compute the mean infection proportion, mean(mean logit-beta), and mean(std dev logit-beta) for the last 10% (last 100,000 steps out of a million).

Haven’t spent enough time to figure out what’s going on here yet …

Mean logit-beta:

Mean proportion infected:

Mean std dev of logit-beta:

Compare high-transmission-equilibrium (red) vs low-transmission equilibrium (blue). Add proportion infected to the plot …

Variable-\(\gamma\) models

Preliminary results for allowing infectious period to mutate: everybody gets infected, all the time! Infectious period (\(1/\gamma\)) just keeps increasing … (becomes very, very long).

I’ve also run the latin hypercube runs for the variable-gamma case, but I’m not sure I want to worry about the results until I’ve understood what’s going on in this case. Related to the details of order/what the effective \(R_0\) etc. is in this particular discret-time model?

To do